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In this paper we study the generation of primordial perturbations in a cosmological setting of 
bigravity during inflation. We consider a model of bigravity which can reproduce the ACDM back¬ 
ground and large scale structure and a simple model of inflation with a single scalar field and a 
quadratic potential. Reheating is implemented with a toy-model in which the energy density of 
the inflaton is entirely dissipated into radiation. We present analytic and numerical results for the 
evolution of primordial perturbations in this cosmological setting. We find that the amplitude of 
tensor perturbations generated during inflation is sufficiently suppressed to avoid the effects of the 
tensor instability discovered in Refs. [1, 2] which develops during the cosmological evolution in the 
physical sector. We argue that from a pure analysis of the tensor perturbations this bigravity model 
is compatible with present observations. However, we derive rather stringent limits on inflation from 
the vector and scalar sectors. 

PACS numbers: 04.50.Kd, ll.10.Ef 


I. INTRODUCTION 

The question whether the graviton may have a mass has attracted considerable attention in the last decade. 
However, constructing a viable theory of massive gravity is a non-trivial problem since the presence of a mass term 
in the gravity action removes diffeomorplrism invariance: the metric has six degrees of freedom (four being absorbed 
by the Bianchi identities), five of these represent the massive graviton while the sixth is usually a ghost, the so-called 
Boulware-Deser ghost [3, 4]. To remove this ghost one needs an additional constraint. This has been achieved with a 
very specific form of the potential for the gravitational field, the dRGT (de Rham, Gabadadze, Tolley) potential [5-7], 
which has been the basis for much of the recent work on this topic (see, e.g., [8-11] and refs, therein). 

In massive gravity theories, the mass term is defined with respect to a fixed reference metric and the possible 
solutions of course strongly depend on this reference metric. Moreover, even when choosing the reference metric to be 
Friedmann, the resulting solutions either do not show the well known cosmological behavior, or they are unstable [12- 
16], see [17] for a review. 

Also from a theoretical point of view, it is rather unsatisfactory to introduce the reference metric as an ‘absolute 
element’, i.e., a non-dynamical held in the theory. For this reason, bimetric (or more general multi-metric) theories, 
with a dynamical reference metric, are more natural [18-20]. Investigations of theoretical aspects of bimetric massive 
gravity can be found in [21-28]. 

Cosmological solutions of bimetric theories can actually fit the expansion history of the accelerating Universe [29- 
33] . Observational tests of several models of bigravity are presented in [34-38] . The cosmology of bigravity in various 
cosmological settings is studied in [39, 40] while in Refs. [41-43] the cosmology of models of bigravity where matter 
is coupled to a combination of the two metrics is investigated. 

Cosmological perturbations in bigravity have been studied in different settings and for different models in [44-47]. 
A more generic study of instabilities in bimetric theories can be found in Refs. [1, 48]. Recently, scalar perturbations 
of these models have been investigated and it has been shown that there exists a class of models of bigravity that 
admit solutions with scalar perturbations free of exponential instabilities at all times, while the other models do 
exhibit exponential instabilities in the scalar sector [40, 49]. 
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The evolution of tensor perturbations in this particular class of models has been first studied in [1] and in more 
detail in [2]. The problem of how cosmological observations are affected by these tensor instabilities and possible 
ways out are discussed in [50]. In [51] a general analysis of the tensor sector in models which are free from known 
instabilities is presented and it is discussed how measurements of the amplitude of primordial gravitational waves can 
be used to constrain them. 

In Ref. [2] it has been found that tensor perturbations are affected by a power-law instability connected with the 
violation of the Higuchi bound [52] in the sector not-coupled to matter. This instability is then transferred to the 
physical sector through the coupling between the two tensor modes. By fine-tuning the amplitude of the unstable 
tensor mode to be highly suppressed with respect to the one of the physical sector at the end of inflation, one can 
achieve that the instability does not show up in the physical metric until today. In [2] this fine-tuning is explicitly 
quantified. The problem of the viability of the model is therefore translated into the question: is the amplitude of 
the uncoupled tensor mode after inflation sufficiently suppressed with respect to the one of the graviton coupled to 
matter? 

In this paper we address this question in detail, i.e., we embed the model of bigravity studied in [2] in an inflationary 
scenario and determine the amplitude and the spectrum of primordial tensor perturbations generated at the end of 
inflation. We find an expression for the ratio between the amplitude of the two tensor modes at the beginning of 
the radiation era as a function of the reheating temperature of the inflationary model. We argue that however, the 
amplitude of the non-physical metric is nearly entirely in the constant mode and the growing mode amplitude is 
much smaller, actually too small to affect the ‘physical metric’, i.e. the metric which is coupled to matter in the way 
discussed in [2]. From this linear analysis, the model is therefore not in conflict with observations. 

We also investigate the vector and the scalar sector. We find that in the model considered, in addition to the 
nearly scale invariant inflaton mode, very large vector perturbations and a very red spectrum of scalar perturbations 
are generated during inflation. The condition for the vector mode not to spoil perturbation theory (i.e. back-react) 
during inflation constrains the scale of inflation substantially. 

While we were working on this, a study of primordial gravitational waves in this model has appeared in [53] in a 
larger context. Our analysis goes beyond the results presented in [53]. We perform an analytical study of primordial 
perturbations in all the sectors. We study in detail, both numerically and analytically, tensor perturbations during 
inflation and reheating. The main results of [53] are confirmed. We also find that tensor perturbations of the second 
metric which are generated during inflation are too small to affect the observable gravitational waves in the physical 
metric. 

The paper is organised as follows. In the next section we present the equations of motion of bimetric gravity for 
cosmological (i.e. homogeneous and isotropic) spacetimes. We then specialise to a model which gives an acceptable 
expansion history. In Sec. Ill we present our model of inflation and reheating and we study its background evolution. 
In Sec. IV we briefly review the perturbation equations of bimetric gravity. The study of tensor perturbations is 
presented in Sec. V and discussed in Sec. VI . In Sec. VII we study the generation of vector perturbations during 
inflation and in Sec. VIII we discuss scalar perturbations. Finally, in Sec. IX we conclude. 

Notation: We set c = h = ^Boltzmann = I- M g = 1/V8ttG = M p ~ 2.4 x 10 18 GeV is the reduced Planck mass. We 
work with the metric signature (—, +, +, +), and we restrict to 4 spacetime dimensions. With • and with ' we indicate 
derivatives with respect to physical time and to conformal time, respectively. The conventions for the bigravity action 
are those of [18-20]. We consider only one of the two metrics coupled to matter, and we restrict to minimal couplings. 
For self-consistency, the Hassan-Rosen bi-metric action and the general equations of motion are given explicitly in 
Appendix A. 


II. COSMOLOGICAL ANSATZ AND BACKGROUND EQUATIONS 

We consider solutions of bigravity where both metrics are spatially isotropic and homogeneous. For simplicity, we 
also assume that both metrics have flat spatial sections, K = 0. Modulo time reparameterizations, the most general 
form for the metrics (in conformal time r) is 

g^dx^dx 1 ' = a 2 (r) (—dr 2 + Sijdx‘dx J ) , (1) 


ffjvdx^dx" = b 2 (r ) (— c 2 (t)cIt 2 + Sijdx l dx J ) . 


(2) 
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Here a and b are the scale factors of the two metrics and c is a lapse function for /. It is convenient to define both 
the conformal Hubble parameter (H) and the standard one (H) for both metrics 



a 


a 

2 


Hf 


b b 2 c ’ 


(3) 


where ' denotes the derivative with respect to the conformal time r. We introduce also the ratio between the two 
scale factors 


r 


: 

a 


(4) 


In the matter sector, we consider the energy-momentum 
equation of state p = wp and 4-velocity u tl . Explicitly, 


tensor of a covariantely conserved perfect fluid with 


= (p + p) +pg^u, 

(5) 

p' = -3 H(p + p), 

(6) 

p = wp. 

(7) 


The general Lagrangian of bimetric gravity and the resulting modified Einstein equations for the metrics g and / are 
presented in the Appendix A. 

The Bianchi constraint in the cosmological ansatz can be written as 


Pg =-3H ( Pg +PG) , 

where we have introduced a ‘gravity fluid’ with density and pressure given by 


( 8 ) 


PG = 


PG = 


87 tG 


2 

TYl* 1 

8nG (& r 3 + 3r 2 + 3/3i r + /3 0 ) , 

(9) 

(/3 3 cr 3 + (2c + l)r 2 + /3i(c + 2 )r + /3 0 ) . 

(10) 


Here the /3j are the parameters of the bigravity potential, see Appendix A. It is easy to show that the Bianchi constraint 
( 8 ) is equivalent to 


m 2 (/?3 r 2 + 2/? 2 r + [5 1 ) {cba! ~ ab') = 0 . 

The equations of motion (the Friedmann equation and the acceleration equation) for the metric g are 

3 H 2 = 8nG {p + PG ) , 


( 11 ) 


( 12 ) 


3 H 2 + = -8irG (p + pc) , 


(13) 


while for the / metric we find the equations of motion 


3„J = ^(A + % + ?A +A 

J \ r 6 r z r 


(14) 




(15) 


where a is a dimensionless parameter, a = Mf/M g (see also Appendix A). Under the rescaling —> a ~ 2 / MJ , and 

/3 n — > a n /3 n , the equations of motion become independent of a [20, 36], which has motivated many works on the 
cosmology of bigravity to simply set a = 1, as we shall do here. Recently, however, it has been argued that this choice 
actually hides the possibility to recover General Relativity (GR) with a cosmological constant in the limit a —> 0, see 
[39] for a detailed discussion. 
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We distinguish two branches of solutions, depending on how the Bianchi constraint (11) is implemented. Either 
there is an algebraic constraint for r, 


/3 3 r 2 + 2/3 2 r + /3i = 0 , 


(16) 


or 


U f =U, rH f = H . 


(17) 


At the background level the first branch with constant r is equivalent to GR with an effective cosmological constant, 
while the second one gives rise to a richer cosmology. We will focus on the second branch in the rest of this work. 
The Bianchi constraint in the second branch can be re-written as 


c = 


r' + rW 
rH 


(18) 


This fixes c as a function of 7~L, r and r'. 

From now on, we will focus on the so-called ‘/ 3 i /34 model’ of bigravity, where all the (3 n parameters but /3i and (3 4 
are set to zero. This model is also called the ‘infinite branch / 3 i /34 model’ or ‘infinite branch bigravity’ in Ref. [49], 
referring to the fact that the initial condition for r has to be chosen in such a way that r evolves from infinity to 
a finite value during the cosmological evolution, in order for the exponential instabilities in the scalar sector not to 
show up. As already mentioned, this model is the only one free of these instabilities. The study of the cosmological 
evolution of this model has been addressed in a series of recent papers [1, 2, 49]. 


III. SCALAR FIELD INFLATION AND REHEATING 
A. General setting 

In this work we focus on the evolution of the model of bigravity during the inflationary period, where the 
dynamics of the universe is dominated by a scalar field <f>, the inflaton, minimally coupled to the physical metric g. 
We consider a simple model of inflation with a single scalar field with mass M$ and quadratic potential. We choose 
the best-fit values pim 2 = 0.48 Hq and /I 4 m 2 = 0.94 Hq obtained in [38] and [49] by fitting measured growth data 
and type la supernovae. 

The Lagrangian density for the inflaton can be written as 

^ = (19) 

The field <fi can in principle interact with other fields such as fermions, gauge bosons, etc., but we assume that this 
interaction can be neglected during inflation and that energy and pressure are dominated by the contribution from 
the inflaton. The energy-momentum tensor of </> is given by 

Tpv =dfj,^d v (j) + g^Ccj, = d v <j> - (^g a0 d a <j) dp(j) + E(</>)^ - (20) 

For the energy density and pressure this yields 

H = ~ T ° = 2^ + 2^ (V0)2 + VW + V(<t>) ~ V{<P) ’ (21) 

and 

p^=^Pt= J f = ^-^ (V0 ) 2 - V(0) V(cj>) ~ -Vtf ). ( 22 ) 

The first approximation in eqs. (21,22) is due to the fact that here we suppose that there exists some (sufficiently 
large) region of space within which we may neglect the spatial derivatives of <j> at some initial time r,;, explicitly 
V<(>(x, ri) < <//(x,Ti). The second approximation is due to the fact that we also suppose that in this region of space 
also the time derivative is much smaller than the potential, <^>(x, t;) <C V 1 ^ 2 These slow-roll conditions are such 
that we have p$ = and p$ + 3 Prf, ~ —2 V(</>) < 0 . 
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At early times in some sufficiently large patch, the Universe is dominated by the potential of a slowly varying (slow 
rolling) scalar field, and hence it is in an inflationary phase. As time goes on, the scalar field starts evolving faster 
and inflation eventually comes to an end when the time derivative of (f> grows to the order of U 1 / 2 . When inflation 
ends, (f) decays rapidly and starts oscillating about its minimum. At the end of inflation, the inflaton oscillates as 

<t> — </>o cos (Mf r) . (23) 


The field <f> is a damped harmonic oscillator with frequency M$. For a harmonic oscillator, when averaging over one 
period we have 


so that 


{V) = 


W 2 ) 

2 a 2 ’ 



and hence (p^) oc a 


(24) 


(25) 


We assume that during these oscillations the coupling of <j> to other degrees of freedom than g^ becomes relevant and 
the inflaton finally decays into a mix of elementary particles which rapidly thermalize. As a simple approximation to 
this complicated and model dependent reheating process, we describe the coupling with the other degrees of freedom 
by means of a dissipative term oc Tcj)' in the equation of motion. In physical time, the equation of motion for the 
inflaton becomes 


p + = -VjM). 


(26) 


During inflation H 3> T and particle production is negligible. When H ~ T, reheating takes place and the inflaton 
energy is rapidly dissipated into other particles which couple to the inflaton. 

We consider a toy-model of reheating in which the energy density of the inflaton is entirely dissipated into radiation. 
In this setting, the total energy momentum tensor has a contribution given by the inflaton and one given by radiation, 


T — qrW _i_ T 1 ( 7 ') 
w -L ■ 

(r) 

Initially X)v = 0. The total energy momentum is covariantly conserved: 

=0 => . 
In our cosmological setting, eq. (28) is equivalent to the following set of equations: 


H + 3 - (P<t> + P<t>) = (P<t> + P4>) 


(27) 


(28) 


(29) 


Pr + 3 (Pr T Vr) — T (p^ “t“ P^) * 


(30) 


It is easy to check that eq. (29) is equivalent to the equation of motion for the inflaton, eq. (26). 


B. Background evolution during inflation 


To study the background evolution during inflation, we consider as a complete set of independent equations the 
two Friedmann equations, (12), (14), and the acceleration equation for the p-metric, (13), the Bianchi constraint in 
the second branch, (17), the equation of state for the inflaton, (22) and the equation of motion for the inflaton, (26). 
Solving the two Friedmann equations together with the acceleration equation for g , we can express r’. H and p r as 
functions of r and p^ 

r r' -3r (1 + w r ) (/34 r 3 - 3/?ir 2 +/3i) + 3r 2 (w r — w^) 87rG?n -2 p^ 

H = n = 2/3 4 r 3 -3/3ir 2 -ft ’ ( ’ 


H 2 



Pi + P 4 r 3 


3 r 


(32) 
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8 n Gp r = m 2 — — 3 m 2 far + m 2 /3 4 r 2 — ^nGp^ , (33) 

r 

where the suffixes V and ‘0’ refer to radiation and to the inflaton, respectively; with p$ = -^2 4>' 2 + ^M 2 (f > 2 and 
ojcj, = Pcf,/P 4 , = —1 + (j)' 2 /(a 2 ptj,). The equation of motion for the inflaton, eq. (26), in conformal time can be written 
as 


<j>" + 2'H<j>' + aT<// + a 2 V,/)((/)) = 0 . (34) 

The two differential equations (31) and (34) are coupled and we solve them together with initial conditions given 
deep in the inflationary epoch. We choose the expectation value of the inflaton at the beginning of inflation to be of 
order 10 M p . Since during inflation the slow-roll condition holds and T -C H, the initial conditions for eq. (34) can 
then be parametrized as 


= 10Mj 


P ’ 


(35) 


P, ,_Ya, \ 

a T * 3H T * 3 H 


(36) 


where we choose for the mass of the inflaton M$ ~ 0.2 eV. 1 2 * Therefore, the state parameter for the inflaton can also 
be written as 


17 


= -1 + 


2 M 2 
3fiim 2 r 2 


(37) 


where in the last equality we have used the fact that during inflation r 1 (as follows from eq. (32)) and p r ~ 0. 
Eq. (32) and r 1 also imply that during inflation 


r(ri) 


3 H(n ) 2 

TO 2 /?4 



(38) 


Once the coupled differential equations (31) and (34) are solved with initial conditions (35), (36) and (38), the 
evolution of the Hubble parameter, of the lapse c and of p r (via eq. (33)) can be derived. 

The results of the numerical integration are shown in Figs. 1, 2 and 3. For the numerical integration, the parameter 
T in eq. (34) has been chosen such that T = H(z re h), where z re h = 5 • 10 1 ' is the reheating redshift. Fig. 1 shows that 
the inflaton starts oscillating at the end of inflation, and that this oscillation is transferred to and c, which starts 
from the value c = 1 during inflation and becomes c = — 1 in radiation domination. The variable r is almost constant 
during inflation (77 ~ 10 33 ) and it starts to decay rapidly in the radiation dominated era. Fig. 2(a) shows that the 
physical Hubble parameter is almost constant during inflation and then starts to decrease. Fig. 3 shows that at the 
end of inflation the energy density of the inflaton is matter-like while the energy density of radiation produced by the 
decaying of the inflaton has the usual evolution with time 2 oc a~ 4 . 


IV. ANALYSIS OF PERTURBATIONS: GAUGE INVARIANT VARIABLES 


We consider perturbations around the Friedmann backgrounds, 


T a h 


g 111* , 


f»v = f, 


flV 


b 2 h 


f ft" ■ 


(39) 


In this section, background quantities are indicated with an overbar. We parametrize the perturbations in as follows: 

-2 A n 


{hg liv) ~ 


C gJ - djBg 

C g i — diB g + diVgj + djV g i + 2 didjE g + 2 SijF g 


(40) 


1 Therefore, from R(n) 2 ~ ^ V(ri) ~ (§^)) 2 it follows that H{n) ~ 1 eV. 

2 We have also checked that the evolution of p r from eq. (33) is equivalent to the one obtained solving the differential equation (30), with 

vanishing initial condition for p r at early times. 
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(a) 



(c) 



(e) 



(b) 



16 18 20 22 


log, 0 (z+l) 

(d) 



logi 0 ( z +D 

(f) 


FIG. 1: We show the ratio between the scale factors r, the lapse of the /-metric c, the equation of state parameter of the 
inflaton and the expectation value of the inflaton as functions of redshift during inflation. We have chosen <j>(zi) = 10 M p 
and H(zi) = 1 eV. The parameter Y in eq. (34) has been chosen such that Y = H(z le h) with z re h = 5 ■ 10 17 . Note how the 
oscillations of <f> lead to strong oscillations of the lapse function c and the equation of state parameter at the end of inflation. 




(a) (b) 


FIG. 2: The physical and the conformal Hubble parameters as functions of redshift, panels 2(a) and 2(b) respectively. 
































FIG. 3: Evolution of the energy density of the inflaton and of radiation (in blue), normalized with respect to the critical 
energy density of the universe today. The yellow curves in panels 3(a) and 3(b) are oc a~ 3 and oc a -4 , respectively. 


with 


(hf 


-2c 2 A f C f j — djB f 

Cfi — diBj fij? + diVfj + djVfi + 2didjEf + 25ijFf 


«<+/= a +,/= = o, 


yih TT 
0 n gJii 


= 0. 


(41) 


(42) 


Spatial indices are raised and lowered using the flat spatial metric <5 tJ . 

In the scalar sector we have 8 fields and 2 gauge freedoms, hence we can form 6 gauge invariant combinations which 
can be chosen as [34, 44] 


* g = A g -nr g A g -r' g , */ = A f + c - 2 (£ - cu f ) r f -c~ 2 r> f , 

$ g = F g -nr g , f = Ff — c _1 Tlf Tf , (43) 

e = E g -E f , B = B f -c 2 B g + (l-c 2 )E ' g , 

where T g j = B g j + E' g In the vector sector we have 4 fields and 1 gauge freedom and hence we can form 3 
gauge-invariant combinations which we choose as follows [44, 47] 


Vgji — C g ,fi Vgji 1 Xi Ggi Cfi . 
The energy-momentum tensor for the perturbed universe is 


(44) 


Tt = Tf* + 8TH . 


(45) 


The perturbations can be divided in perfect-fluid and non-perfect-fluid ones, with 5+5 dof (degrees of freedom). The 
perfect fluid dof in ST/f are those which keep Tjf in the perfect fluid form: 


Tt = (P + p)u^u u +p5^ . 


(46) 


We suppose here that the perturbations are only of this type. Thus, they are given by the density perturbation, the 
pressure perturbation and the velocity perturbation. Explicitly: 

p = p + Sp , p = p + Sp , u l = u l + 5u l = Su l = -Vi. (47) 

a 

The Su° is not an independent dof, it is fixed by the normalisation u tl u u g^ v = — 1. 

We can now write the perturbed Einstein equations for the two metrics. In the following we will use the Fourier 
transform of perturbations with respect to x\ the corresponding 3-momentum will be k’ and k 2 = To keep the 

notation simple, the Fourier transform will be denoted by the same symbol as the original function. 


V. TENSOR PERTURBATIONS 


Tensor perturbations of a given k-mode are composed of two independent helicity modes, 


lTT 


= fr + 4 +2) + h~ 


(- 2 ) 

ij 


(48) 
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where + and — denote the two helicity-2 modes of the gravitational wave. For an orthonormal system (k, e^\e* 2 )) 
we have 


= ± 4 (1) ± 

V 


ie 


(2) 


and 


J+2) 


— e + e + 

— • 


0 (- 2 ) 


In what follows we assume parity invariant perturbations, 


(/ l + (k)(/ l + (k / ))*} = k 3 (h-(k)(h-(k')y) = s{k - k , )27r 2 p h (k ), 


(49) 


and ( h + h ) = 0. And we shall consider just one mode, say /ij = hfGf and /i+ = h g G g , where Gf and G g are 
uncorrelated Gaussian random variables with vanishing mean and variance (G g j(k)G g j(k')) = d(k — k')27r 2 , so 
that h g j is the square root of the power spectrum. All the following is also valid for the modes h~ * which are not 
correlated with j in the parity symmetric situation which we consider. 

With a perfect fluid source term, i.e. no anisotropic stress, in the first order modified Einstein equation, we obtain 
the following tensor perturbation equations for our bimetric cosmology [2]. 

h” + 2 Hh' g + k 2 h g + m 2 a 2 r /3i (h g - hf) = 0 , (50) 


2 

h'f + c 2 k 2 hf — m 2 Pi^— (h g — hf) = 0 . (51) 

In Ref. [2] we have solved these coupled differential equations in the radiation era and have found that hf has a 
growing mode, hf oc r 3 on large scales which via the coupling enhances also the mode h g of the physical metric. Here 
we solve these equations numerically and analytically in the inflationary regime, where sensible approximations can 
be introduced to simplify the system. 


2 (H + - ) - - 

r 1 c 


A. Analytical results during inflation 

Deep in the inflationary epoch, the potential V(</)) is very flat and the inflaton is slowly rolling. Since pf, — —pf,, 
it is legitimate to model this period as a de Sitter phase with constant Hubble parameter H = Hj ~ const, (where 
the suffix T hereafter stands for inflation). From eq. (32) it follows that during inflation r = 77 = const., with 
r 2 ~ 3H 2 /(m 2 (34 ) ~ 3 (Hj/H 0 ) 2 and eq. (18) gives for the lapse function of the /-metric c ~ const ~ 1. With the 
parametrization T~L = —1/r and a = — 1/(tHj) (note that with this choice r < 0 during inflation), and m 2 /3ia 2 ~ 
( H 0 / Hj) 2 t~ 2 . eqs. (50) and (51) can be approximated in a de Sitter universe as 

K ~ 2 -K + k2 b 9 + (I 7 ) 4 ~ h s) = 0 - (52) 


h"j ~\h' f + k 2 h f ^ ^ = (53) 

These equations can be solved exactly in terms of oscillating and decaying modes. 

We want to choose as initial conditions the quantum vacuum of the graviton degree of freedom. For tensor 
perturbations the canonically normalised variables (recall that Mf = M g = M p ) are given by 

(Qg)ij = dj Q g = M p a (hfj r ) g , (Q f ).. = e« Qf = M p b (hJ g T ) f . (54) 

Equations (52) and (53) in terms of these new variables and recalling that b = r a become 

oi '+(* 2 - Q ‘ + {w) h (°»“ \ Qt ) = 0 ■ (55) 

0/ + ( t2 -^) Q >-{%) (5 = (56) 
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Since during inflation 77 ~ Hi/Hq 1, eqs. (55) and (56) can be approximated as 


Q'a 


k r 2 ) Qg + 


Ho 

H! 


2 Qg - 0, 


(57) 


«/-(§) 


(58) 


For sub-horizon scales, |/cr| 1, eqs. (57) and (58) reduce to two copies of the same equation for a harmonic oscillator 

with frequency k. The quantum vacuum solutions are 


Q g = ^ / _ exp (—i hr ), Qf = ex P * ^ r ) > f° r |fcr| ^ 1 • 


(59) 


We want to solve eqs. (57) and (58) with initial conditions (59). These equations can be decoupled introducing the 
new variable Q + = Qf + ( 777 ) Qg 


q 9 + (|n ^q 9 = 0 , 


(60) 


Q'\ 


k ‘ -yi ) 0 + = °- 


(61) 


Eqs. (60) and (61) can be solved in terms of Bessel functions. Requiring that the asymptotic behavior (59) is recovered 
for |Air| 1, we find the following solutions for the canonically normalized variables Q g and Qf 


Qg — 


p r I hr T . . pv jkr . . 

V 2 S2k J i^ {kr)+ 'P \l2k r h^^- 


(62) 


Qf = 


p2k 


l + w, 


1 fcr 1 6 


—i kr 


-(^1 Q„- 


(63) 


To simplify the notation, we have introduced the tiny constant 7 = Hq/Hj , 3 In the limit 7 —t 0 we have Q g = Qf = 

Q+- 

The canonically normalized variables Q g and Qf are connected to the power spectrum by 


p hg {k) = k 3 I hi? hi jTT \ = 2 ■ fc3 IQff 


*7 9 9 


a 2 


(64) 


P hf{k )=k 3 \hP^ f hJ TT \ = ^ 


kpQf? 

a 2 MS ’ 


(65) 


where the factor of 2 is due to the two hclicity modes in each tensor sector. Hence from eqs. (62) and (63) we can 
find the solutions of eqs. (52) and (53) making use of the relations 


h n = - 

9 a M. r 


k 3/2 Q q , 


h f = —- \y- k 3 t 2 Qf . 

1 rj aM p 1 


( 66 ) 


3 Given that for \kr\ 1, the behavior of the Bessel functions is J 1 v / 9 _ 47 (A:r) —> — yj cos (kr) and Yi v / 9 _ 47 (Air) — > — yj sin (kr), 

the asymptotic behavior of (62) and (63) for \kr\ 1 is exactly of the type (59). 
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B. Numerical results during inflation and reheating 


The asymptotic behavior of the solutions (66) for \kr\ 1 during inflation is given by 


h 


a 


k 

V2 M p 


Hi t e 


—i kr 


hf = - 


k 

y/2 M p 


Hit c - ikT 
ri 


for \kr\ 1. 


(67) 


These functions and their first derivatives can be evaluated at r = n to And the initial conditions for the numerical 
evolution of the full tensor perturbation equations, (50) and (51). The results of the numerical integration are shown 
in Figs. 4 and 5, for four different fc-modes. 

Independently of the mode k, the agreement between the analytical solutions for h g and hf, obtained from (62) 
and (63), and the numerical one is reasonably good in the regime in which the slow-roll condition holds and the 
background is well approximated by de Sitter. As expected, both modes h g and hf oscillate in the redshift range in 
which k > H(z), see eq. (67). Furthermore, the mode hf develops an instability at the end of the inflationary period, 
where it oscillates with increasing amplitude. This is due to the fact that the damping term in eq. (51) becomes an 
anti-damping term at the end of the inflationary stage. Indeed, using eq. (18), eq. (51) can be written as 


h" f 


2 cH- 


h'f + c 2 k 2 hf — m 2 /?i- (h g — hf) = 0 . 


( 68 ) 


Since c = 1 at the beginning of the inflationary era whereas c = —1 in the radiation era, the term in square bracket 
changes sign when inflation ends, going from 2 H to —2 H, i.e. from a positive damping term to a negative anti-damping 
term. 

From eqs. (62) and (63), it follows that on super horizon scales the power spectra at the end of inflation are scale 
invariant and given by 


— r / Ph f (z, k ), |fcr|<l. (69) 

This result has also been derived in [53]. Ph g hence has the same behaviour as the standard (i.e., GR) tensor power 
spectrum, whereas Ph f is suppressed with respect to Ph g by a huge factor r 2 . The numerical results for the power 
spectra at the end of inflation are shown in Fig. 6. In the analytical result (69) slow roll corrections are neglected 
since in this context we are mainly interested in orders of magnitude and not in very precise results. The power 
spectra shown in Fig. 6, however, are numerically calculated with the full model. 


VI. DISCUSSION: IS THE MODEL STILL VIABLE? 

In [2], the cosmological evolution of tensor perturbations in the model of bigravity has been addressed and the 
condition needed for the instability not to show-up until present times has been quantified in terms of a fine-tuning of 
the amplitude of the two tensor modes after inflation. We have found that during the radiation dominated era tensor 
fluctuations of the / metric can grow like a 3 on super horizon scales. Furthermore, if they become larger than those 
of the g metric, h/(to) > h g (to), they can influence the latter via the coupling term and show up e.g. in the CMB. 

Considering naively, as from eq. (69), hf/h g (z e ) = 7'J 1 ~ Hq/Hi at the end of inflation and that hf grows by an 
amount (T eq /T e ) 3 ~ ((1 + z eq )/(l + z e )) 3 , the condition for the instability not to show-up hf(t 0 ) < h g (to) [2] implies 
a very low bound on the value of Hi and therefore on the scale of inflation. Here the index e denotes the value of a 
quantity at the end of inflation. 

However, this naive argument is misleading since at the end of inflation, when \kr e \ <C 1 for the relevant modes, 

hf ~ r ^w ( x + i fcTe i 2 )' ( 70 ) 

Here, r e ~ — (a e f//) _1 ~ —(1 + z e )/Hi is the conformal time at the end of inflation. At the beginning of the radiation 
era the constant mode of hf turns into the constant mode and only the very severely suppressed decaying mode, 
oc |fcr e | 2 turns into the growing mode. Therefore, considering only this decaying mode for hf at the end of inflation 
(when we put our initial conditions on the amplitudes of the tensor modes), our bound will be reduced by a factor 
| kT e | 2 ~ ((1 + z e )k/Hi) 2 . This bound is valid for perturbations on very large scales with wave number k ~ H 0 . 
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(a) 



(c) 



(e) 

k=2 10 13 Ho 



(g) 



(b) 



(d) 



(f) 

k=2 10 13 H 0 



(h) 


FIG. 4: Evolution of tensor perturbations of the metrics g and / as functions of redshift. The numerical solution (blue) is 
plotted together with the analytical one (in red) valid in the inflationary era, i.e., in the regime in which the hypothesis of 
slow-rolling holds. We have chosen k = 10 1() Ho , k = 10 11 Ho , k = 10 12 Ho and k = 2-10 13 Ho in the panels 4(a)-4(b), 4(c)-4(d), 
4(e)-4(f) and 4(g)-4(h), respectively. The spectrum for the <?-mode is rescaled with a factor r/ ~ 10 33 while the spectrum for 
the /-mode is rescaled with a factor rj. Note that in our model inflation ends roughly at log 10 (l + z) ~ 19.0 while radiation 
domination is established at log 10 (l + z) ~ 17.5. 
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(a) 

k=10 n Ho 



(c) 



(e) 



(g) 


k=10 10 H 0 



(b) 


k=10 n Ho 



(d) 



(f) 


k=2 10 13 Ho 



16 18 20 22 


l°gi 0 (z+ 1 ) 

(h) 


FIG. 5: We show the same plots as in Fig. 4 but in log scale for more detailed apprehension of the decay and growth of 
perturbations. One sees that the analytical solution during inflation is out of phase with the numerical solution. The reason 
is that the space-time background is somewhat different in the two cases, in fact, for our analytical solution the background is 
pure de Sitter, whereas for the numerical solution we have taken into account the full evolution of the background. 
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(a) 



(b) 


FIG. 6: Power spectra at the end of inflation, z — 10 19 . The power spectrum for the hf mode has been rescaled by r?, with 
77 ~ 10‘ i9 for easier comparison with the spectrum of the physical mode h g . 


Inserting this value for k , and replacing ((1 + z eq )/( 1 + z e )) 3 


(Hi/H eq ) 3 / 2 by (H T /H 0 ) 3 / 2 for simplicity 4 , we obtain 


¥i «»>« 






1/2 

< 1 . 


(71) 


For the last ~ sign we used that r 7 1 ~ Hq/Hi. Therefore, no meaningful bound is obtained. Also for the smallest 
scales, with wave number k — Hi/( 1 + z e ), there is no significant bound since inside the horizon the growth of hf is 
only linear in a so that the factor ( T eq /T e ) 3 has to be replaced by ( T eq /T e ) and the same result is obtained. 

One may finally ask how model independent is the fact that the dominant constant mode at the end of inflation 
only transits into the constant mode of the radiation era. On super horizon scales and for negligible couplings the 
mode equations (50) and (51) are always of the form 


h" + a(r)h' = 0, 

so they always have a constant mode. The general solution of (72) in fact is simply 


h(r) = A\ J dr exp 


a(T n )dT" 


+ A 2 . 


(72) 


(73) 


The first mode is decaying or growing, depending on the sign of a while the second mode is always constant. Hence 
even if there is a relatively long reheating phase, there will also be a constant mode during this phase. 

Therefore, from a pure analysis of tensor perturbations, the model cannot be ruled out. 

Note also that our inflationary model is typically above the strong coupling scale, A 3 , of the theory [54], For 
T re h ~ 10 MeV we have Hi ~ T? eh /Mp < 1CT 22 GeV ~ A 3 = (m 2 Mp ) 1 / 3 ~ 2 x 10 ~ 22 GeV. Therefore, if the 
reheating temperature is above lOMeV, bigravity becomes strongly coupled and it is not granted that cosmological 
perturbation theory still applies. However we underline that, since the strong coupling scale is derived in a Minkowski 
background, it is not entirely clear whether this represents an upper limit on the Hubble scale or just on the energy 
of the perturbations on a given background. For the sake of the argument, our approach here is to simply analyse the 
classical bigravity Lagrangian in a Friedmann Universe with small perturbations. 


VII. VECTOR PERTURBATIONS 

Vector perturbations of a given k-mode can be decomposed as 

V i = V (1) ef ) +V (2) ef ) , (74) 

where the two orthonormal vectors eand are defined in Sec. V. In what follows we shall consider just one 
mode, say since the situation is perfectly symmetric for the other mode and suppress the superscript, so that 


The last replacement is an approximation which corresponds to consider radiation domination lasting until present time. 
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Vi = e* V. If the background is pure de Sitter, in the vector sector only the mode Vi = V g i — Vfi propagates [44]. The 
action for this vector mode in Fourier space can then be written as 5 


P / dTdJt F+aWft OKI 2 - (fc 2 + aWA) IV. I 2 ) , 


Sv = 


where e.g. |Vi| 2 = V*V‘. The canonically normalized variable in this case is defined as 


Q i =e l Q = M p a 2 k J ' ™ ^ 

y k. + a z r m. p i 


V,;. 


After integration by parts, the action (75) for the variable Q can be written as 

Sv = J dr d 3 k ^ [(Q') 2 - C(k,r) Q 2 

where, in order to simplify the notation, we have introduced 


(75) 


(76) 


(77) 


C(k , r) = k 2 + Pi m 2 r a 2 — 2H 2 — %' 


Pi m 2 r a 2 + k 2 


-3 w 


Pi m 2 r a 2 + k 2 


(78) 


Using that pirn 2 ~ Hq and that in pure de Sitter Universe with Hubble constant Hj, a = —1 /(Hit) and 77 ~ Hi/H 0l 
the previous expression can be approximated as 


C(A,t)=^+[^ - 2 - 


Ho \ 1 / (fcr) 2 \ 1 / ( kr ) 2 \ 3 2 2 6 

^ " Uo/ff; + N 2 j ^ " \H 0 /H I + (kT) 2 J r 2 ~ ^ ~ ^ : 


(79) 


where in the last equality we have used that Hq/Hj <C 1 and we have assumed that also Piin 2 a 2 r <C k 2 holds, 
or equivalently Hq/Hj <C (kr) 2 . This second inequality is valid for the following reason: during the radiation era 
ra 2 =const.~ -\/3f2 ra d (This can be derived from the two Friedmann equations, see [2].). Since this quantity is growing 
during inflation we have r/a 2 < -\/3U rad . Therefore Pim 2 ra 2 ~ H^rja 2 ~ ( H 0 /Hj)t~ 2 < 77g\/3H rad < k 2 for all 
values k>H 0 which are observable. 

The equation of motion derived from the action (77) with C(k , r) ~ k 2 — is then 6 


Q"+lk 2 --^) Q = 0. 


(80) 


For | kr\ 1, Eq. (80) reduces to a harmonic oscillator equation with frequency k, and has the vacuum solution 

1 


Q = 


—ikr 


\/2k 


\kr\ » 1. 


(81) 


Eq. (80) can be solved exactly. Asking that the asymptotic behavior (81) is recovered for \kr\ 1 the solution is 
given by 


Q = 


ikT v(2) 


\/2k 


K ’{kr ), 


(82) 


where lip is the spherical Hankel function of the second kind of order l, see [55]. Substituting eq. (82) in eq. (76), 
the evolution of the physical variable V can be written in terms of Q. Using again that 77 ~ Hj/Hq, pim 2 a 2 r <C k 2 
and that Pirn 2 ~ Hq, we obtain from eq. (76) 


1 1 1 


I Hr 


Vi ~ M v a 2 H t V H 0 Qi ' 


(83) 


5 This action is equivalent to the action (63) in ref. [47]. 

6 One can also verify that the exact equation of motion for Q which can be derived from the action (77) with the exact expression for 
C(k , t) coincides with eq. (7.9) in [44], once written in terms of the original variable V. 
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For superhorizon modes \kr\ <C 1, using the asymptotic behaviour of the spherical Hankel function, h^\x) ~ —3/a; 3 
for \x\ -C 1, we find 


V, ~ 


3 Hj 

75 k 



d . 


Therefore, for super-Hubble scales, the vector power spectrum can be written as 


Pv(k,r) 



(HA* 

\Mj H 0 ’ 


\kr\ <C 1, 


(84) 


(85) 


where the multiplication by a factor 2 in the last expression is due to the two vector modes and the powers of k are 
introduced to make the power spectrum dimensionless. Note the large enhancement by a factor Hi I Ho with respect 
to the standard tensor spectrum which is of order ( Hj/M p ) 2 . 

In order for linear perturbation theory to remain viable, one has to request at least that P\> < 1, which means that 
our inflation model must be such that Hj < 10 -2 GeV (where we have used the fact that Ho ~ 10~ 42 GeV and that 
M p ~ 2.4 • 10 18 GeV). This requires a rather low scale of inflation which is however acceptable. 

Asking that vector perturbation should not be larger than scalar perturbations after inflation, P\> < 10 -9 [56] 
would require an inflationary Hubble scale of Hj < 10 -5 GeV which corresponds to a reheating temperature of 
T reh < (ff/Mp) 1 / 2 ~ 10 6 GeV. This inflationary scale is not excluded, see [57, 58] but rather low. However, we have 
not studied the evolution of vector perturbation during the radiation era. If they decay, as in GR, this second limit 
is not relevant. It only applies if vector perturbations in bigravity stay constant during the radiation dominated 
Universe. 


VIII. SCALAR PERTURBATIONS 

Let us now turn to scalar perturbations. For this we assume that, like for the other degrees of freedom, the difference 
to GR during inflation mainly comes from the existence of additional modes, but that the GR-modes are not strongly 
affected, since the coupling between the GR-modes and the additional modes is suppressed by Ho/Hj. We therefore 
assume that the inflaton mode leads to a nearly scale invariant spectrum like in GR, and we only study the additional 
helicity-0 mode of the massive graviton. For this we work in a pure de Sitter background and neglect the slow roll and 
the inflaton perturbation. In this situation the hclicity-0 mode of the massive graviton is the only dynamical scalar 
degree of freedom. It is given by a linear combination of the two Bardeen potentials [44], 

$ = - 2 rj$ f . ( 86 ) 


Its evolution is governed by the equation 

<F" + 2 


2 k 4 


9a 2 U 2 m\ + k 4 - 18 H 4 


- 1 


l* 


4 (fc 6 - 3k 4 H 2 ) 
9a 2 U 2 m\ + k 4 - 18H 4 


+ 3a 2 m| — k z — 6H~ 


= 0, 


(87) 


where 



m 2 Pi ~ Hq . 


( 88 ) 


Using the same approximations as for the vector mode, eq. (87) can be approximated on sub-Hubble scales by 


d*" + 277$' + fc 2 $ = 0 , |fcr|»l. 


(89) 


Analogously to tensors, we quantize the scalar perturbations in order to find the initial conditions. The canonical 
variable is given by <j> = M p a<&. In terms of this variable, eq. (89) reduces to a harmonic oscillator equation with 
vacuum solution 4>(t) = e~ lkT /V2k. It follows that 


ff/ e“ 

Mp 77 k 


\kr\ 1. 


(90) 
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On super-Hubble scales eq. (87) can be approximated by 

$" - 2H$' - = 0 , |fcr| <C 1, 


(91) 


with general solution 


d> = Cit +, |fcr|<l, (92) 

where Ci and C 2 are integration constants, which can be fixed by matching the sub-Hubble solution and its derivative 
with the one in the super-Hubble regime. 

Note that the mode oc C 2 manifests an instability on super Hubble scales since |r| is decreasing during inflation. 
This is the manifestation of the fact that also during inflation the Higuchi bound is violated in the scalar sector. 
Indeed, for the scalar sector (helicity-0 mode) the Higuchi bound of the /3 i/ 34 model is given by [1, 33] 

3r + - 2/3 4 r 2 > 0, 



which is violated for r > 1.02. In our treatment, neglecting couplings of the scalar mode to other modes, the instability 
coming from this violation only shows up as a growth of perturbations on super-Hubble scales which leads to a red 
spectrum as we now show. Note also that the growth is exponential in physical time since r -2 oc a 2 oc exp(2 Hit). 
Working out the matching conditions explicitly we obtain 


$ = —e 


H, 1 
M p ^2 k 



i 1 \ 

T + 3 t 2 fc 3 ) ’ 


\kr\ -C 1, 


(93) 


The power spectrum for scalar perturbations on super-horizon scales can then be expressed as 

P*(T,k) = fc 3 |$| 2 ~ {jf) (t) ocfc ” 4 ’ 1^1 <!■ ( 94 ) 

This very red power spectrum is strongly enhanced on large scales, |fer| <C 1. Comparing it to the standard 
inflationary scalar power spectrum which is of the order of [56] 

where e < 1 denotes the slow roll parameter, one must conclude that this mode, if it transits to the radiation era 
completely spoils the observed large scale structure. However, for scalar perturbations the matching from inflation to 
the radiation era has to be studied carefully, it can even lead to a change in the power spectrum as found, e.g., for 
the inflationary magnetic mode studied in Ref. [59]. For this reason, we shall not draw strong conclusions from this 
result. 

Nevertheless, we request that P$(z,k) < 1 for perturbation theory to remain valid during inflation, so that we 
can neglect back-reaction of the perturbation to the cosmic evolution. At the end of inflation we have r/a 2 nd — 
(rn 2 )rad — V3H r so that Hend = 17"end| — 1 = H/a en d ~ (HiHq) 1 / 2 (3H r ) x ! 4 . Inserting fc ~ Ho in ('^ond/fc) 4 we obtain 
(Hend/Ho) 4 ~ 3f l r (Hj/H 0 ) 2 which leads to 


('T'end) H 0 ) 


3H r ( Hj \ 
18^ \M P H 0 ) 


The condition P$(r en d,Ho) < 1 then becomes 


Hi < 


18 MlH% 


3n r 


1/4 


10’ 11 GeV, 


Vj /A ~ (H/Mp) 1 / 2 < 10 4 GeV. 


(95) 


(96) 


Also this is indeed a rather low inflation scale. Requesting < 10 9 would reduce it by another two orders of 
magnitude. 
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IX. CONCLUSIONS 

In this paper we have studied the generation of perturbations during inflation in a bimetric theory of gravity. We 
have analysed the evolution of the two tensor modes and we have found that both acquire a scale invariant spectrum 
with hf = hg/ri, where the ratio 77 = (b/a)\j — Hj/Hq 3> 1 is nearly constant during inflation. In addition to this 
constant tensor mode, the /-metric sector has also a decaying mode, which is the one that turns into the growing 
mode during the subsequent radiation era. Nevertheless, this mode is so severely suppressed that it does not lead 
to a significant amplification of the physical tensor mode as discussed in Ref. [2]. Therefore, looking at the tensor 
sector alone, the Pi Pa model of bimetric gravity cannot be ruled out despite the fact that the Higuchi bound of the 
tensor sector of the /-metric is violated. Note that this Higuchi bound on a Friedmann background does not lead to 
an exponential instability but only to power law growth of fluctuations. For this reason, the detailed analysis of the 
initial conditions from inflation carried-out in this work was needed to decide whether the model is ruled out or not. 

We have also analysed the vector (helicity 1) and scalar sectors. Also vector perturbations are generated during 
inflation leading to a scale invariant vector spectrum with an amplitude which is boosted by a factor 77 with respect 
to the tensor spectrum. Requiring vector fluctuations to remain perturbative gives an upper limit to the scale of 
inflation, Hj < 10~ 2 GeV which translates to an inflationary energy scale < 10 8 GeV. In principle, even if the 
scale of inflation is lower, these vector mode will source tensor modes in the /-metric at second order which then can 
feed into the growing mode. Anyway, a detailed calculation of this is beyond the scope of the present work. 

In the scalar sector we have not discussed the inflaton perturbations, assuming that they are not modified due to 
the very weak coupling of bigravity during inflation. However, in a bigravity theory we have the helicity-0 mode of 
the graviton as a second scalar mode. As the Higuchi bound in the scalar sector is violated this mode is growing 
on super Hubble scales during inflation. We have computed its spectrum at the end of inflation and have found 
that it is very red, oc k~ 4 . Requesting that these fluctuations remain perturbative also on the largest scales k ~ Hq 
gives stringent constraints on the scale of inflation, Hi < 10 -11 GeV, which translates to an inflationary energy scale 
V/ /4 <10 4 GeV. 

We conclude that the Pi Pa model studied in this paper cannot be ruled out from the analysis of the tensor sector 
alone. Nevertheless, it is strongly constraint due to the large vector perturbations which are generated during inflation 
and due to the very red spectrum of scalar perturbations. However, considering that all the problematic scales of 
inflation have Hi > 10“ 11 GeV, which is far higher than the strong coupling scale, A 3 , these results have to be taken 
with a grain of salt. Still, we recall that this strong coupling limit is derived in a Minkowski background and it is not 
clear that it should invalidate quantum perturbations on classical Friedmann solutions with a Hubble scale which is 
larger than A 3 . 

All other models of bigravity where matter only couples to one of the metrics (the g metric in this work) suffer from 
exponential instabilities of scalar perturbations on a FLRW background. This makes these models less attractive as 
candidate solutions to the dark energy problem. Due to the breakdown of linearity one has to work out the theory 
at higher orders and hope to cure the instabilities, possibly through the Vainshtein mechanism [60]. A possible 
way out is to push the gradient instability to very early times, rendering it unobservable. This can be achieved by 
lowering the value of the Planck mass of the metric which does not couple to matter [39]. In addition, there remain a 
multitude of massive (bigravity) models whose cosmology deserves further investigation, e.g., where matter, or even 
different matter sectors, can couple to both metrics [23, 25, 61, 62]. Alternatively one could also consider non-FLRW 
backgrounds or change the status of the parameters of the theory, e.g. by promoting the Pi coefficients to functions 
of the helicity-0 mode [28] , or of some independent scalar field. 
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Appendix A: Hassan-Rosen bigravity model: general aspects 

The conventions used for the bigravity action are those of [20]. Only one of the two metrics is coupled with matter, 
and we restrict to minimal couplings. The action is given by 


S = — / d 4 X\f^g 


M 2 

R(g)-2m 2 V(gJ)) + C m (g^) 


f . , _ M 2 

- / d 4 x\f-~]-J-R(f ), 


(Al) 


where g is the physical metric (the one coupled to matter), / is the second metric, and M g = 1/\/8 ttG = M p and Mf 
are the respective Planck masses with dimensionless ratio a = Mf/M g . We assume the matter fields $ to be coupled 
to g only. The potential is given in terms of the tensor field X = \Jg~ 4 f : 


V(g,f) = '£ l p n e n (X). 


n—0 


where the coefficients (3 n are constants and the polynomials e n (X) are 


e o — I; e i — [X]) 

e 2 = i([X] 2 -[X 2 ]), 

e 3 = p:([X] 3 — 3[X][X 2 ] + 2[X 3 
o 


1 


e 4 = — ([X] 4 -6[X] 2 [X 2 ]+8[X][X 3 ]+3[X 2 ] 2 ^6[X 4 ]) = detX. 
The square bracket [• • • ] denotes the trace. The equations of motions for g pl , and / M „ are 


(A2) 

(A3) 

(A4) 

(A5) 

(A6) 


R + ^ £(-)" & [sma Y* )u (n/PP) + 9ux (\/PP) 


n —0 




T 

^ Liv j 


(A7) 


2 . 


f^R+ T 2 ^ E(-) n Pi-n [Ux F ( *)„ (PPP) + fuX (PPP) 


= 0 , 


n —0 


where the overbar indicates f pv curvature. The definition of the (X) matrices is as follows: 

V( 0) (X)=I, T ( i) (X) = X - I [X] , 

Y(2) (X) = X 2 - X [X] + li 
y (3) (X) = x 3 - x 2 [X] + lx ([X] 2 - [x 2 ]) - ([X] 3 - 3 [X] [x 2 ] + 2 [x 3 ]) . 


(AS) 

(A9) 

(A10) 

(All) 


As a consequence of the Bianchi identity and of the covariant conservation of T, lv , we find the following Bianchi 
constraints (for each one of the two metrics) 




n—0 


9n\ Yf n)v (\/pP) + 9v a Y^ (y/PP) 


= 0, 


(A12) 


V M 


6 

£(-)"&-„ [pa Yfc )v (PPP) +f»x Ev (PPP 


n—0 


= 0 , 


(A13) 


where the overbar indicates covariant derivatives with respect to the / metric. Both these constraints follow from the 
invariance of the interaction term under the diagonal subgroup of the general coordinate transformations of the two 
metrics. They are equivalent and in this work we focus on the first one. 
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